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1 Introduction 

Spatial optical solitons are known to originate from the nonlinearity-induced 
diffraction suppression and beam self-trapping in a bulk dielectric medium . 
Since, generally speaking, the effect of diffraction is strong, considerable op- 
tical nonlinearities are required in order to compensate for the diffraction- 
induced beam spreading. As a result, in the vicinity of a self-trapped beam 
the refractive index experiences large deviations from a Kerr-type dependence 
and the models of generalized nonlinearities, that describe beam self-trapping 
phenomena and spatial solitons, become nonintegrable. 

For solitary waves of nonintegrable nonlinear models, linear stability is one 
of the crucial issues, since only stable (or weakly unstable) self-trapped beams 
can be observed in experiment. The study of soliton stability in nonlinear op- 
tical materials has a long history, also associated with the study of nonlinear 
waves in other media such as plasmas and fluids. Stability of one-parameter 
solitary waves has already been well understood for both fundamental (single- 
hump and nodeless) solitons |^J^,^,^ and solitons with nodes and multiple 
humps 1^,0. The pioneering results of N. Vakhitov and A. Kolokolov 
known these days as the Vakhitov- Kolokolov stability criterion, found their 
rigorous justification in a general mathematical theory developed by M. Gril- 
lakis et al. Although the corresponding stability and instability theorems 
for the scalar nonlinear Schrodinger (NLS) models formally extend to the 
case of multi-parameter solitons most of the examples analyzed so far 
correspond to solitary waves with a single parameter. Only very recently a 
systematic analysis of more general cases has been carried out, in connection 
with the study of three- wave parametric solitons in quadratic (or x'^') optical 
media 0. 

Recent progress in the study of soliton instabilities can be associated with 
the application of a multi-scale asymptotic bifurcation theory developed for 
weakly unstable stationary nonlinear localized waves. In the framework of this 
theory, the unstable eigenvalue of the associated linear spectral problem is 
treated as a small parameter of the asymptotic expansions [ p^ , and the re- 
sulting equation for the slowly varying soliton propagation constant may also 
account for weakly nonlinear effects that describe the long-time evolution of 
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linearly unstable solitary waves, and an effective saturation of the soliton in- 
stability due to higher-order nonlinear effects. In the case of multi-parameter 
solitary waves, a simplified version of this asymptotic method applied near 
a marginal stability point (or, in general, a surface) is reduced to finding 
certain determinants constructed from the derivatives of the system invari- 



ants ]ll| , |12| , |l3 14 . However, the validity of this bifurcation theory has no 
rigorous proof, and it can only be used to estimate the domains of the soli- 
ton stability and instability. Since more general oscillatory instabilities may 
also occur |7|, [l5| , |l6|jl^ jl^ , numerical simulations are often required in order 
to verify the predictions of the asymptotic theory (see, e.g., Ref. as an 
example) . 

Recently, a general matrix criterion for the stability and instability of 
multi- component solitary waves was derived ||l9|| for a system of N incoher- 
ently coupled NLS equations. In this general approach, the soliton stability 
is studied as a constrained variational problem reduced to finite-dimensional 
linear algebra. Unstable eigenvalues of the linear stability problem for multi- 
component solitary waves were shown to be connected with negative eigenval- 
ues of the Hessian matrix constructed for the energetic surface of A^— component 
spatially localized stationary solutions. 

The analysis and, correspondingly, stability criteria obtained in Ref. 
can be extended, at least in principle, to other types of solitary waves, such 
as incoherent solitons in non-Kerr media, parametric solitary waves in x^^-* 
optical media, etc. In all such cases, the results on stability and instability 
of solitons can be readily obtained with a rigorous generalization of some 
of the previously known results of the multi-scale asymptotic theory. How- 
ever, in each of those cases some additional analysis is required, in order to 
clarify whether the results of the asymptotic theory completely define the 
stability properties of multi-component solitary waves. Beyond the validity 
of the multi-scale analysis, oscillatory instabilities may occur, and appropri- 
ate studies should rely solely on the numerical analysis of the corresponding 
eigenvalue problems. 

In this Chapter, we present a brief overview of the basic concepts of 
the soliton stability theory and discuss some characteristic examples of the 
instability-induced soliton dynamics, in application to spatial optical solitons 
described by the NLS-type models and their generalizations. First of all, we 
demonstrate a crucial role played in the stability theory by the soliton internal 
modes that present an important characteristic of solitary waves in nonin- 
tegrable nonlinear models. In particular, we study an example model of the 
NLS equation with higher-order nonlinearity, in order to show that the soliton 
internal mode gives birth to the soliton instability. Near the marginal stabil- 
ity point, the soliton stability can be analyzed by a multi-scale asymptotic 
technique when the instability growth rate is treated as a small perturbation 
parameter. We also discuss some results of the rigorous linear stability analy- 
sis of fundamental solitary waves and nonlinear impurity modes. More recent 
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studies of higher-order sohtons revealed that the so-caUed multi-hump vector 
solitary waves may become stable in some nonlinear models, and we discuss 
the stability of (l+l)-dimensional two- and three- hump composite solitons 
created by incoherent interaction of two optical beams in a photorefractive 
nonlinear medium. The final part of this Chapter is devoted to the stability 
of solitons in higher dimensions and, in particular, it discusses very recent 
results on the symmetry-breaking instability of a (2-|-l)-dimensional vortex- 
mode composite soliton and the formation of a rotating dipole-like structure 
associated with a robust radially asymmetric dipole-mode vector soliton, a 
new composite object resembling a "molecule of light" . 



2 Linear eigenvalue problem 

To discuss the stability properties of spatial optical solitons, we consider 
the nonintegrable dimensionless generalized NLS equation that describes the 
(l-l-l)-dimensional beam self- focusing in a waveguide geometry, 

i^ + g + ^(/;.)^ = 0, (1) 

where "0(0;, z) is the dimensionless complex envelope of the electric field, x is 
the transverse spatial coordinate, z is the propagation distance, / = \ip{x, z)\'^ 
is the beam intensity, and the real function J-^{I] x) characterises both linear 
and nonlinear properties of a dielectric medium, for which we assume that 
J^(0; ±oo) = 0. Stationary spatially localized solutions of the model (|]) have 
the standard form, ip{x, z) — <P{x; f3)e^^^ , where (3 is the soliton propagation 
constant (/3 > 0) and real function (!>{x\ (3) vanishes for \x\ oo. An impor- 
tant conserved quantity of the soliton in the model (|l|) is its power defined 
as 

\^p{x,z)\^dx ^ / <P^{x](3)dx. (2) 

-oo J —oo 

To find the linear stability conditions, we consider the evolution of a small- 
amplitude perturbation of the soliton presenting the solution in the form 

i^ix, z) = |<Z>(a;; /?) + [v{x) - w{x)] e'^' + [v*{x) + w*(a;)]e-*^*^} e''^^(3) 

where the star stands for a complex conjugation, and obtain the linear eigen- 
value problem for v{x) and w{x), 

Lqw ~ Xv, Liv — Xw, 

d^ (4) 
dx"^ 



where Uo = x) and Ux = x) + 2I[dT{I; x)/dl]. 
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A stationary solution of the model is stable if all the eigenmodes of 
the corresponding linear problem do not have exponentially growing am- 
plitudes, i.e. Im(A) = 0. It can be demonstrated that the continuum part of 
the linear spectrum of the problem (Q) consists of two symmetric branches 
corresponding to real eigenvalues with the absolute values |A| > /3, and there- 
fore only discrete eigenstates are responsible for the stability properties of 
localized waves. Then, the corresponding eigenmode solutions fall into the 
following categories: 

• internal modes with real eigenvalues describe periodic oscillations; 

• instability modes correspond to purely imaginary eigenvalues; 

• oscillatory instabilities can occur when the eigenvalues are complex. 

In what follows, we present several approaches that allow us to study 
analytically and numerically the structure of the discrete spectrum in order 
to determine the linear stability properties of solitary waves. We also analyze 
the nonlinear evolution of unstable solitons. 



3 Soliton internal modes and stability 

Since the soliton instabilities always occur in nonintegrable models, it is in- 
teresting to know what kind of distinct features of the solitary waves in 
nonintegrable models might be responsible for their instabilities. It is com- 
monly believed that solitary waves of nonintegrable nonlinear models differ 
from solitons of integrable models only in the character of the soliton interac- 
tions: unlike "proper" solitons, interaction of solitary waves is accompanied 
by radiation |2^. However, the soliton instabilities are associated with non- 
trivial effects of different nature that are generic for localized waves of nearly 
integrable and nonintegrable models. In particular, a small perturbation to 
an integrable model may create an internal mode of a solitary wave pl| . This 
effect is beyond a regular perturbation theory, because solitons of integrable 
models do not possess internal modes. But in nonintegrable models such 
modes may introduce qualitatively new features into the system dynamics 
and, in particular, lead to the appearance of the soliton instabilities. 

To demonstrate that internal modes are generic for nonintegrable models, 
we consider a weakly-perturbed cubic NLS equation with the nonlinear term, 

T{I:x)^I + efiI), (5) 

where /(/) describes a deviation from the Kerr nonlinear response, and e is a 
small parameter. Then, the stationary solution can be expressed asymptot- 
ically as ^(a;;/3) = ^o(a^) + -I- 0(6^), where <Po{x) = y/2j3sech {^/J3x) 
is the soliton of the cubic NLS equation, and 'Pi{x) is a localized correction 
derived from Eqs. and Neglecting the second-order corrections, we 
find the results for the effective potentials of the linearized eigenvalue prob- 
lem (|), t/o = <?o + ^Uo and Ui = 3^>g + eUi, where Uq = f{<Pl) + A^q^i 
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Fig. 1. Schematic presentation of the origin of the bifurcation- induced soliton 
instabilities: (a) spectrum of the integrable cubic NLS model, (b) bifurcation of 
the soliton internal mode, (c) collision of the internal mode with the neutral mode 
resulting in the Vakhitov-Kolokolov-type soliton instability 



and Ui = f{<Pl) + 2(Plf'{<Pl) + 12<?o^i, and the prime denotes diflFerentiation 
with respect to the argument. 

The linear eigenvalue problem (^) and (||) can be solved exactly at e = 
(see, e.g., Ref. |^). Its discrete spectrum contains only the degenerated 
eigenvalue at the origin, A = 0, corresponding to the so-called soliton neutral 
mode [see Fig. |l|(a)]. We can show that a small perturbation can lead to the 
creation of an internal mode (that corresponds to two symmetric discrete 
eigenvalues), which bifurcates from the continuous spectrum band, as shown 
in Fig. 0(b). To be definite, we consider the upper branch of the spectrum 
and suppose that the cut-off frequencies Am = ±/3 are not shifted by the 
perturbation. Then, the internal mode frequency can be presented in the 
form A = /3 — e^K^, where k is defined by the following result pl| ]: 

\K\^jsign{e) f {v{x,P)UiV{x;f3) + W{x,f3)UoW{x;P)\dx. (6) 

Here {V{x; f]),W{x; P)} are the eigenfunctions of the cubic NLS equation 
calculated at the edge of the continuous spectrum, V{x] (3) = 1—2 sech^(v^a;) 
and W{x]l3) = 1. A soliton internal mode appears if the right-hand side of 
Eq. (^ is positive. 

As an important example, we consider the case of the NLS equation (|l|),(|) 
perturbed by a higher-order power-law nonlinear term, 

f{I)^^l'- (7) 
The first-order correction to the soliton profile can be found in the form 
V2/3'"^/2 [2cosh(2V;5x) -t- cosh{A^/^x)\ 



<?i(x) 



3cosh^(v^a;) 
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Fig. 2. Soliton instability in the model (§),(H) and (^), presented through (a) the 
power dependence -P(/9) and (b) the evolution of the discrete eigenvalue of the 
problem (3) that defines the soliton internal mode (solid) and an instability mode 
(dashed). Solitons for /3 > /3cr, i.e. for dP/dj3 < [dashed curves in (a) and (b)] 
are linearly unstable. Dotted lines show the asymptotic dependences calculated 
analytically 



With the help of Eq. (||) , it is easy to show that for e > a perturbed NLS 
sohton possesses an internal mode that can be found analyticahy near the 
continuum spectrum edge, 



A = /3 



(8) 



For high intensities, the additional nonhnear term (0) is no longer small, 
and the soliton solutions, together with the associated hnear spectrum, should 
be calculated numerically. Power dependence P{(3) calculated with the help 
of Eq. (||) for the soliton of the model (|l|),(||), and (Q) is presented in Fig. ^(a), 
and it is matched with the discrete eigenvalue of the linearized problem 
shown in Fig. ^(b). First of all, we notice that the asymptotic theory provides 
accurate results for small intensity solitons, i.e. for /3 < 0.1 (shown by the 



Stability of Spatial Optical Solitons 



7 



(o) (b) 



z 




Fig. 3. Evolution of a perturbed unstable NLS soliton in the model (§),(H),(0) for 
e = 1, and /3 = 2, in the case of (a) increased power (collapse) and (b) decreased 
power (switching to a low-amplitude stable state). The initial power was changed 
by 1% compared to the exact soliton solution 



dotted curves in the insets). Secondly, the slope of the power dependence 
changes its sign at the point (3 = /?cr, where the soliton internal mode vanishes 
colliding with the soliton neutral mode, as depicted in Figs. 0(b,c). At that 
point, the soliton stability changes due to the appearance of a pair of unstable 
{purely imaginary) eigenvalues [dashed curve in Fig. ^(b)]. In Sec. ^ below, 
we prove rigorously a link between the soliton stability and the slope of the 
dependence P(/3). 

For [3 ~ Pen the instability-induced dynamics of an unstable soliton can 
be described by approximate equations for the soliton parameters derived 
by the multi-scale asymptotic technique (see Sec. ^ below) but, in general, 
we should perform numerical simulations in order to study the evolution of 
linearly unstable solitons. In Figs. ^(a,b), we show two different types of the 
instability-induced soliton evolution in our model. In the first case, a small 
perturbation that effectively increases the soliton power results in an un- 
bounded growth of the soliton amplitude and subsequent beam collapse [see 
Fig. H(a)]. In the second case, a small decrease of the soliton power leads 
to a switching of a soliton of an unstable [dashed. Fig. ||(a)] branch to a 
stable [solid. Fig. ||(a)] one, as is shown in Fig. ||(b). This latter scenario be- 
comes possible because, in the model under consideration, all small-amplitude 
solitons are stable. However, if the small-amplitude solitons are unstable, 
the soliton beam does not converge to a stable state but, instead, diffracts. 
Thus, in the NLS-type nonlinear models there exist three distinct types of the 
instability-induced solitons dynamics [ p^ . 
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4 Stability criterion for fundamental solitons 

Direct investigation of the eigenvalue problem (Q) is a complicated task which, 
in general, does not yield a complete analytical solution. However, for a class 
of "fundamental" solitary waves (i.e. solitons with no nodes), the analysis can 
be greatly simplified. First, we reduce the system (^) to a single equation: 

LqLiv = X^v, (9) 

for which the stability condition requires all eigenvalues to be positive. It is 
straightforward to show that Lq = L'^L~, where = zLd/dx +<P~^ {dip / dx) , 
and thus, instead of (|^), one can consider an auxiliary eigenvalue problem psf , 

L-LiL+i = X^i, (10) 

which reduces to Eq. (||) after the substitution v = L+w. Since the operator 
L~LiL+ is Hermitian, all eigenvalues A^ of Eqs. (||) and ( |lo| ) are real, and 
in this case oscillatory instabilities do not occur. 

Properties of the operators Lj {j — 0, 1) are well-studied in the literature, 
in particular, as a characteristic example of the spectral theory of the second- 
order differential operators (see, e.g., Ref. Q). For our problem, we use 
two general mathematical results about the spectrum of the linear eigenvalue 
problem LjipiP = aI^V""*: 

• the eigenvalues can be ordered as A„_[_]^ > A„ , where n > defines the 

(i) 

number of zeros in the corresponding eigenfunction (fn , 

• for a "deeper" potential well, Uj{x) > Uj{x), the corresponding set of the 
eigenvalues is shifted "down", i.e. An < An . 

Let us first discuss the properties of the operator Lq , for which the soliton 
neutral mode is an eigenstate, i.e. Lq(1>{x; (3) = 0. As we have assumed earlier, 
^(x; /3) > is the ground state solution with no nodes and, therefore, > 
Aq""* = for n > 0. This means that the operator Lq is positive definite on the 
subspace of the functions orthogonal to <?(a;;/3), which allows to use several 
general theorems (^||J^,0,^ in order to link the soliton stability properties to 
the number of negative eigenvalues of the operator Li . Specifically, 

• the soliton instability appears if there are two (or more) negative eigen- 
values, i.e. A^^'' < 0; 

• the solitons are always stable provided the operator Li is positively def- 
inite; 

• in an intermediate case, the soliton stability depends on the slope of the 
power dependence P{P), according to the Vakhitov-Kolokolov stability 
criterion Q|, i.e. the soliton is stable if dP/d/S > 0, and it is unstable, 
otherwise. 
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Thus, to distinguish between these cases, it is sufficient to determine the signs 
of the zeroth and the first eigenvalues of the operator Li. 

To demonstrate the validity of the Vakhitov-Kolokolov criterion, we follow 
the standard procedure ||]. First, wc note that, for the fundamental solitons 
with no nodes, both the direct, Lq, and inverse, Lq^, operators exist, and 
they are positively definite for any function orthogonal to <l^(x] which can 
only be a neutral eigenmode of Eq. (0), and therefore can be ignored. Thus, 
by applying the inverse operator Lq to Eq. (|^), we obtain another linear 
problem with the same spectrum: 

Liv = X'^Lq^v, (11) 

where v(x) now satisfies the orthogonality condition: 

r + oo 

{v\<P)= v*{x)^x](3)dx ^0. (12) 



Second, we multiply both sides of Eq. ( |ll[ ) by the function v*{x), integrate 
over X, and obtain the following result. 



2 {v\Liv) 



{v\Lq v) 

Because the denominator in Eq. (|l^) is positively definite for v satisfying 
Eq. (|l|, the sign of this ratio depends only on the numerator. The instability 
will appear if there exists an eigenvalue < (so that A is imaginary), and 
this is possible only if 

min {v\Liv) < 0, (14) 



where we normalize the function v{x) to make the expression in Eq. ( |14[ ) 
finite, 

{v\v) = 1. (15) 

In order to find the minimum in Eq. (|l^) under the constraints given by 
Eqs. and (|l|), we use the method of Lagrange multipliers, and look for 
a minimum of the following functional: 

C = {v\Liv) - V {v\v) - ^ , 

where real parameters v and [i are unknown. With no lack of generality, we 
assume that /i > 0, as otherwise the sign of the function v(x) can be inverted. 
The extrema point of the functional L can be then found from the condition 
5C/5v — 0, where 5 denotes the variational derivative. As a result, we obtain 
the following equation: 



Liv = vv + ii<P, 



(16) 
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where the values of v and /i should be chosen in such a way that the condi- 
tions dl^ ) and (^5|) are satisfied. Then, it follows from above that {v\L\v) — 
V {v\v) and, according to Eq. (p^, the stationary state is unstable if and only 
if there exists a solution with < 0. 

Operator L\ has a full set of orthogonal eigenfunctions ip'n' [|2^, i-e. 

^'n^\^'m'^ = if n 7^ TO. The spectrum of L\ consists of discrete (a1^^ < 
/3) and continuous (A^^-* > /3) parts, and the norms of the eigenmodes, 

^^n' \H>^n''^ 1 are scaled to unity and a delta- function, respectively. Then, we 
can decompose the function vix) in the following way: 

where the sum goes only over the eigenvalues of the discrete spectrum of ii. 
Coefficients in Eq. ( p7| ) can be found as £>„ = (^Lpn''\v'^. Function ${x; f3) can 

be decomposed in a similar way, with the coefficients C„ = (^(pn^^ . Then, 
Eq. (|l6|) can be reduced to: 

M C„/(a1'^ - u), if Cn y6 and /i > 0, 
Dn = <( 1, if Cn = 0, /i = 0, and = xi^\ (18) 
0, otherwise. 

Note that we should not consider degenerate solutions v = 0, and therefore 
yU = is only possible if v = Xn'' and C„ = 0. In order to find the La- 
grange multiplier ly, we substitute Eqs. (17) and ( p^ ) into the orthogonality 
condition ([T^), and obtain the following equation for the parameter v. 

r+oo 

Q{u) ^ («|<P) - + / CnDldX^^'> = 0. (19) 

As has been mentioned above, instability appears if there exists a root v < 
0, and thus we should determine the sign of minimal v solving Eq. (|l9|). 
Because the lowest-order modes of the operators Lq and Li, <P{x] (3) and i^^^ 
respectively, do not contain zeros, the coefficient Cq ^ 0. Then, from the 
structure of Eq. (|l^ ) it follows that Q(v < Aq^^) > 0, and thus the solutions 
are only possible for v > Aq^"*, meaning that if Aq^'' > the stationary state 
${x;(3) is stable. 

We notice that the function Q{v) is monotonic in the interval (— oo, +oo) 
for /X > and Aq"'^'' < v < X^n \ where n > 1 corresponds to the smallest 
eigenvalue with C„ 7^ 0. Then, it immediately follows that instability appears 
if A^^^' < and Ci 7^ 0. On the other hand, if Ci — 0, the corresponding 
eigenmode ip"i^ satisfies Eq. (^6|) and the constraints (p^,(p^) with v = A^^'' 
and fjL = Q and, therefore, an instability is always present if X^-^^ < 0. 
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The last possible scenario is when Aq^'' < and A^^'' > 0, so that the modes 
with V = An^'' do not give rise to instabihty, and we search for the solutions 
with ^ > 0. Then, because a1^^ > A^^^ > 0, the sign of the solution v is 
determined by the value of Q(0). Indeed, if (5(0) > 0, the function Q(y) 
vanishes at some < which indicates instability, and vise versa. From 
Eqs. ® and (|l|), it follows that Q{v 0) = (ij^V^I^)- To calculate this 
value, we differentiate the equality Lq(1> = with respect to the propagation 
constant and obtain 

d<P 

i^.^ = -«^, (20) 

that finally gives: sign[(3(0)] = —sign{dP/df3). These results provide a proof 
of the stability conditions outlined above. 

An important case is the soliton stability in a homogeneous medium, when 
does not depend on x. Then, a fundamental soliton has a symmetric 
profile with a single maximum, and d<P/dx is the first-order neutral mode of 
the operator Li, i.e. A^^-* = 0. Therefore, in such a case the soliton stability 
follows directly from the slope of the dependence P{P)- 

Linear stability discussed above should be compared with a more gen- 
eral Lyapunov stability theorem which states that in a conservative system 
a stable solution (in the Lyapunov sense) corresponds to an extrema point 
of an invariant such as the system Hamiltonian, provided it is bounded from 
below (or above). For the NLS equation, this means that a soliton solution 
is a stationary point of the Hamiltonian H for a fixed power P, and it is 
found from the variational problem 5{H + j3P) = 0. In order to prove the 
Lyapunov stability, one needs to demonstrate that, for a class of spatially 
localized solutions, the soliton realizes a minimum of the Hamiltonian when 
P is fixed. This fact can be shown rigorously for a homogeneous medium with 
cubic nonlinearity Q|, i.e. for !F{I\x) = /, and it follows from the integral 
inequality, 

H > H, + {P^/^ ~ P}^^), (21) 

where the subscript "s" defines the integral values calculated for the NLS 
soliton. Condition ( pl| ) proves the soliton stability for both small and finite- 
amplitude perturbations, and a similar relation can be derived for generalized 
nonlinearity, being consistent with the Vakhitov-Kolokolov criterion (see de- 
tails in Ref. [|). 



5 Marginal stability point: Asymptotic analysis 

As follows from the linear stability analysis, the solitons in a homogeneous 
medium are unstable when the slope of the power dependence is negative, 
i.e. for dP/d(3 < 0. Near the marginal stability point (3 — /3cr defined by 
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the condition {dP/d(3)p=p^^ = 0, where the instabiUty growth rate is smaU, 
we can derive a general analytical asymptotie model which describes not only 
Hnear instabihties but also the nonlinear long-term evolution of unstable soli- 
tons. Such an approach is based on a nontrivial modification of the soliton 
perturbation theory Indeed, the standard soliton perturbation theory is 
usually applied to analyse the soliton dynamics under the action of exter- 
nal perturbations. Here we should deal with a qualitatively different physical 
problem when an unstable bright soliton evolves under the action of its "own" 
perturbations. As a result of the development of such an instability, the soli- 
ton propagation constant varies slowly along the propagation direction, i.e. 
[3 = (3{z). As the instability growth rate is small near the threshold (3 = (3cn 
we can assume that the profile of the perturbed soliton is slowly varying with 
z and the soliton evolves almost adiabatically (i.e., it remains self-similar). 
Therefore, we can develop an asymptotic theory representing the solution to 
the original model (||) in the form ip — <j){x; (3; Z) exp[z/3o2 + «e /g^ l3{Z')dZ'], 
where (3 — (3o + e^n{Z), Z = ez, and e <C 1. Here the constant value f3o is 
chosen in the vicinity of the critical point (3cr- Then, using the asymptotic 
multi-scale expansion in the form 4> — "^{x; (3) + e^(j)3{x; P; Z) + O(e^), we 
obtain the following equation for the soliton propagation constant /3 (details 
can be found in Rcfs. ^d^), 



M{(3, 



d^f2 1 dP 



dZ^ e2 df3 



2~d^ 



= 0, (22) 

/3=/3cr 



where P{/3) and M{(3) are calculated through the stationary soliton solution, 
and 



M(/3) = 



dx > 0. 



A remarkable result which follows from this asymptotic analysis is the fol- 
lowing. In the generalized NLS equation (|l|) the dynamics of solitary waves 
near the marginal stability point /3 = /3cr can be described by a simple 
collective-coordinate model (^2|) which is equivalent to the equation of motion 
of an effective {inertial and conservative) particle of the mass M{(3ct) with 
the coordinate fi moving under the action of a potential force proportional 
to the difference Pq - P(/3), where Pq = P{Po)- 

First two terms in Eq. ( p^ ) give the result of the linear stability theory, 
according to which the soliton is linearly unstable provided dP/dP < 0. Non- 
linear term in Eq. (|2^) allows to describe not only linear but also long-term 
nonlinear dynamics of an unstable soliton and, moreover, to identify qualita- 
tively different scenarios of the instability-induced soliton dynamics near the 
marginal stabihty point, as discussed in Refs. [jl0[|25|]. 



6 Nonlinear guided waves and impurity modes 
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6.1 Model and general remarks 

In this section, we demonstrate how the basic concepts and results of the 
soliton stability theory can be employed to study the stability of nonlin- 
ear guided waves in inhomogeneous media. In particular, we consider the 
model with a real function J-{I] x) that describes both nonlinear and 
inhomogeneous properties of an optical medium. This kind of problems has 
a number of important physical applications ranging from the nonlinear dy- 
namics of solids to the theory of nonlinear photonic crystals and waveguide 
arrays in optics. In application to the theory of guided electromagnetic waves, 
the model we discuss below describes a special case of a stratified (or layered) 
dielectric medium for which nonlinear guided waves and their stability have 



been analyzed during the last 20 years |p6|j27 |. 

Following the original study presented in Ref. p8[ | , we consider a simplified 
case when the inhomogeneity is localized in a small region, and it is created by 
a thin layer embedded into a nonlinear medium. Then, if the corresponding 
wavelength is much larger than the layer thickness, in the continuum-limit 
approximation the layer response can be described by a delta-function and, 
therefore, we can write 

T{I-x)^F{I)+5{x)G{I), (23) 

where the functions F{I) and G{I) characterise the properties of a bulk 
medium and the layer, respectively. The model (0),(^) describes, in partic- 
ular, a special case of a more general problem of the existence and stability 
of nonlinear guided waves in a stratified medium with a Kerr or non-Kerr 
nonhnear response (see, e.g., some examples in Refs. [p6|). 

When a thin layer is introduced into the system, the translational in- 
variance of the model is broken at the point of the layer location, i.e. at 
X = 0. Nonlinear modes of such an inhomogeneous model should be found as 
spatially localized solutions of the following equation. 



l3^+'^+T{^'^;x)^^Q. (24) 



We notice that for the problem at hand the profiles of localized modes should 
satisfy a homogeneous equation on both sides of the layer and, therefore, they 
can be constructed by using the solutions of Eq. (HJ) with !F{I;x) — F{I). 
Due to a translational symmetry, such a localized solution can be presented 
in the form 'Po{x — xq), where Xq is an arbitrary shift. We also note that the 
profile is symmetric, <l>o{x) = <l>o{—x), it does not contain zeros, ^o{x) > 0, 
and it has a single hump since d<l>o/d\x\ < at a; ^ 0. Therefore, a general 
solution of inhomogeneous Eq. ( p^ ) for a spatially localized normal mode can 
be presented in the following form, 

^.(a:) = (j°J"-^°)',"^0' (25) 
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where xq and s are defined through the layer parameters. From the continuity 
condition at a; = 0, it follows that s = ±1. Thus, the parameter s defines 
the type of symmetry of the localized mode, i.e. the mode is symmetric, for 
s = — 1, and asymmetric, for s = +1. Parameter xq is defined by the following 
matching condition, 

(l-s)^(xo) = G(/o)<Po(xo), (26) 

where Iq = <?g(a;o). We note that the power of the asymmetric modes is the 
same as that for solitons in a homogeneous medium, Pq = 'pQ{x)dx, 
while the power of symmetric modes (hereafter denoted as P) is smaller than 
Pq if Xq < 0, and larger otherwise. 

As follows from Eq. (|25|), the function ^(x) is sign definite. Then, in order 
to apply the general stability results derived for the fundamental solitons of 
the NLS-type equations in Sec. ^, we should study the spectral properties 
of the linear operator Li. Due to the presence of a layer, the number of 
negative eigenvalues of the operator Li depends on the mode parameter /3. 
As has been demonstrated in Ref. |^ , if an eigenvalue passes though zero and 
changes its sign, the so-called critical point appears in the power dependence 
P(/3) calculated for a localized mode. One type of such critical points is 
a bifurcation, where the mode belongs simultaneously to the two families 
of localized solutions, symmetric and asymmetric ones. Thus, the stability 
properties of these distinct types of nonlinear localized modes can be expected 
to be different, and we discuss these two cases separately. 



6.2 Stability of symmetric and asymmetric modes 



Symmetric modes. In this case it can be verified that, similar to the solitary 
waves of homogeneous models (see Sec. ^ above), the operator Li has a 
neutral mode d(I>/d\x\ with zero eigenvalue, but only for a special value of 
the layer response: 



\ dx 



dx^ 



(27) 



where Gi = G(/o) + 2Io[dG/dI]i=ig. For a:o > the function ^{x) has two 
humps, and its derivative d(l>/d\x\ corresponds to the second-order mode with 
A2^-' — 0. On the other hand, for xq < 0, the neutral mode describes a 
ground-state eigenfunction with Aq^'' = 0. From the spectral theorem |Q, it 
follows that this eigenvalue increases for Gi < Gf, and decreases otherwise. 
Additionally, we note that due to the symmetry of the potentials in the 
linear eigenvalue problem, Uj{x) = Uj{~x), the amplitude of the first-order 
eigenmode vanishes at the layer location (x = 0), and thus its eigenvalue A^^^ 



does not depend on Gi. However, the inequality Aq"^ < Xq''' < X'jf' should be 



i(2) 
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fulfilled for any Gi, and we immediately come to the conclusion that A-j^ < 0, 
if xq > 0, and A^^"* > 0, if xq < 0. It is also straightforward to verify that 
for cco = we should have G(/o) = 0, and the derivative d(l>{x)/dx, which 
has a single zero at cc = 0, is a neutral mode of the linear operator Li, i.e. 
A^^'' = 0. In Table ^ we summarize the properties of the operator Li and 
the corresponding general conditions for the stability of symmetric nonlinear 
localized modes. 



Table 1. Conditions for stability of symmetric modes 



Conditions 




Li spectrum 


Stability 


xo > 




< A(i' < 


unstable 


xo < 0, Gi 


< Gf 


< aW < aW 


stable 


xo < 0, Gi 


> GY, dP/dp > 


A<'^ < < aW 


stable 


xo < 0, Gi 


> G^^ dP/d(3 < 


A*'^ < < aW 


unstable 



The next important step is to connect the spectral characteristics of the 
linear operator Li (and thus the stability of nonlinear localized modes) with 
the power dependence -P(/3)- First, we notice that at the bifurcation point 
defined as P(/3) = Po(P), the parameter Xq passes through zero and the 
eigenvalue A^^-* changes its sign (see Table |l|). In particular, for P{(3) > Pa{P) 
such a mode should have a double-hump profile and, therefore, it is always 
unstable. In contrast to these unstable double-hump modes, a family of asym- 
metric modes emerges at the bifurcation point. 

Asymmetric modes. Although asymmetric modes have the profiles of soli- 
tons, the corresponding spectrum of the operator Li can be different. Only 
in the special case Gi = 0, there exists the first-order neutral mode with 
A^^^ = 0, i.e. Li{d$/dx) = 0. 

There can exist several families of asymmetric modes, each of them charac- 
terized by the constant intensity at the layer /q, so that sign(Gi) ~ sign[dG/dI]i^ig 
at the point where the layer response vanishes, G(/o) = 0. If Gi is negative, 
the stability of the localized mode is the same as the stability of a soliton 
in a bulk medium. On the other hand, the modes corresponding to positive 
Gi are unstable, and they can be attracted to or repelled from the nonlinear 
layer. Performing the analysis similar to that in Ref. |^, we find that for the 
asymmetric modes Gi > 0, if after the bifurcation point the branch of the 
symmetric modes is lower than that of the asymmetric modes, P(/3) < Po{/3). 
In the opposite case, we have Gi < 0. In Table 2 we summarize the stability 
conditions for asymmetric nonlinear localized modes. 
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Table 2. Conditions for stability of asymmetric modes 



Conditions 


Li spectrum 


Stability 


Gi > 


< < 


unstable 


Gi < 0, dPo/dp > 


aW < < aW 


stable 


Gi < 0, dPo/df3 < 


aw < < aW 


unstable 



6.3 Example: Competing power-law nonlinearities 

To discuss a characteristic example of the stability theory briefly outlined 
above, we consider the power-law nonlinear response of both the bulk medium 
and the layer, 

F{I)=r, G{I)^a + hP, (28) 

where the nonlinearities are characterized by the positive powers a and 7, 
and the coefflcients a and h account for the linear and nonlinear properties 
of the layer, respectively. 

We consider the case when the layer itself possesses repulsive nonlinearity 
(i.e., 6 < 0). Due to a competition of localizing and de- localizing nonlinear- 
ities, the power dependence of the family of localized symmetric modes can 
become multi-valued, i.e. two or three states may correspond to the same 
value of the mode propagation constant (see Fig. Moreover, such states 
can exist below and above the linear cut-off. 

The branch of the asymmetric solutions emerges at the bifurcation point, 
/3f,, where the power of the symmetric localized modes, P(/3f,), coincides with 
the power of solitons in a homogeneous medium, Po(/3f)) (open circle in Fig. 
Stability of asymmetric modes for 6 < is the same as that of the solitons 
in a bulk medium, i.e. the modes are stable for cr < 2 and unstable oth- 
erwise. On the other hand, after the bifurcation point, i.e. for (3 > P^, the 
symmetric localized mode becomes two-humped and, therefore, it is unstable 
to asymmetric perturbations. In Fig. ^, we illustrate the evolution of a two- 
hump symmetric mode when the power of a bulk nonlinearity is below the 
instability threshold (cr < 2). In this case, the symmetry-breaking instability 
is observed, and the mode is transferred to one side of the layer. Such an 
excitation of a stable asymmetric state is accompanied by a slowly decaying 
quasi-periodic oscillations of the mode amplitude. 

7 Multi-component solitary waves 

7.1 General theory: N coupled NLS equations 

As has been already mentioned in the Introduction, the Vakhitov-Kolokolov 
stability criterion is valid for one-parameter solitary waves, such as scalar 
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Fig. 4. Power vs. propagation constant diagram for (t = 1,7 = 2, a = 2, and 
6 = — 1. Solid and dashed parts of the main curves correspond to stable and un- 
stable localized modes, respectively. Open circle — the bifurcation point from the 
symmetric (S) to asymmetric (AS) modes. Thin dashed — the soliton power Po(/3). 
The modes corresponding to the marked points (a),(b), and (c) are shown below 




Fig. 5. Evolution of a slightly 
perturbed symmetric local- 
ized mode shown in Fig. ^c). 
An unstable two-hump mode 
eventually transforms into a 
stable asymmetric mode 
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non-Kerr solitons, two- wave parametric solitary waves in a quadratic medium, 
etc. As we demonstrated in Sec. ^, the Vakhitov-Kolokolov criterion can also 
describe, in some cases, the stability of nonlinear guided waves in a scalar 
model. However, when solitary waves described by nonintegrable nonlinear 
models possess many parameters, the stability analysis becomes more in- 
volved, and in many cases a direct numerical study of the corresponding 
eigenvalue problem is required. Nevertheless, in some cases it is possible to 
extend the stability analysis to multi-parameter solitary waves. Below, we fol- 
low the recent original work by Pelinovsky and Kivshar ||l9|] and demonstrate 
how the multi-parameter generalization of the Vakhitov-Kolokolov stability 
criterion can be derived for a system of N coupled NLS equations. Such a 
theory includes, as a limiting case, the bifurcation analysis near the marginal 
stability curve. 

We consider a system of N incoherently coupled NLS equations, 

+ dnVlijn Inmli^rnf'^ = 0, (29) 

where stands for Laplacian in the D— dimensional space x = (xi, xu), 
and all the coefficients d„ are assumed to be positive. When one of the vari- 
ables of the vector x stands for time, Eqs. ( |29| ) describe the spatio-temporal 
dynamics of self-focused and self-modulated light in the form of the so-called 
light bullets. 

Provided the symmetry conditions 7„„i — jmn are satisfied, the sys- 
tem (^9|) conserves the Hamiltonian, 

l-oo f N N N \ 

/ dxK]d„|VxV'n|'-^EE^»™l^"l'l^'"l' ' 

\n=l ^„=lm=l / 

the individual mode powers. P„ — \ ! IV'nP'^x, and the total field momentum. 
Localized solutions of Eqs. (|29| ) for the multi-component fundamental solitary 
waves are defined in a standard form as V'n = ^^n(x)e*''"^, where <?„(x) are 
real functions with no nodes, and /3„ are the positive propagation constants. 
Such localized solutions can be also obtained as stationary points of the 
Lyapunov functional, 

N 

A[,l,]=H['^]+Y.Pr.Pnm. (30) 
n=l 

and, therefore, the first variation of A vanishes, whereas the second variation 
6'^A[ip] defines the stability properties: the negative directions of the second 
variation correspond to unstable eigenvalues in the soliton stability problem. 

The stability problem is defined by minimizing the second variation of 
the Lyapunov functional 

/•oo 

5'^A= / dx[(u|Liu) + (w|Low)], (31) 
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where u(x) and w(x) are perturbations of the multi-component soUtary wave 
taken in the form ■0 = $(x) + [u + iw] (x)e'*'^, and the scalar product is 
defined as (f |g) = ^n=i fn9n- The matrix Sturm-Liouville operator Lo has 
a diagonal form with the elements 

N 
m=l 

and the matrix operator Li has the elements 



N 

{Ll)nn = -dnV^ + ^ XI "T"™^™ " '^Jnn'^l, 

m—1 



at the diagonal, and {Li)nm — —'^Inm^n'l'm, off the diagonal. Similar to 
the stability problem for one-component solitary waves discussed in Sees. |^ 
and ^, the matrix operators Lq and Li define the linear eigenvalue problem 
for stability of multi-component solitary waves, 

Liu = —Aw, Lqw — Au. (32) 



Both the linear problem (|3^) and minimization problem (|3l| ) should satisfy 
a set of N additional constraints, 

/oo 
dx('^„e„|u) = 0, (33) 
-oo 

where e„ is the rfi^ unit vector, which correspond to the conservation of the 
individual powers P„ under the action of a vector perturbation (u, w). 

First of all, we recall that the one-parameter solitary waves with no nodes 
(A^ = 1) &re stable in the framework of the constrained variational problem 
(^-(^ provided the energetic surface As{f3) = A[^] is concave up, i.e. 

Under this condition, the linear eigenvalue problem ( |3^ ) has no unstable 
eigenvalues, i.e. those with a positive real part A. Otherwise, the second vari- 
ation (^l|), constrained by the N conditions (^3[), has a single negative direc- 
tion that corresponds to a single positive eigenvalue A in the linear eigenvalue 
problem ( ^ ) (see Refs. |^,|| and discussions above). The stability criterion 
for scalar (or one-component) NLS solitons holds when the self-adjoint oper- 
ator Li has a single negative eigenvalue, i.e. when the second variation (^, 



without the constraint (33) imposed, has a single negative direction. If the 
latter condition is not satisfied, as it happens for solitary waves with nodes, 
the fundamental criterion for the soliton instability can be extended only 
for special cases, while more generic mechanisms of oscillatory instabilities. 
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associated with complex eigenvalues of the linear eigenvalue problem, may 
appear beyond the prediction of the fundamental criterion. 

In some particular cases, the soliton stability analysis can be extended to 
multi-component solitary waves. In particular, this is possible for a system 
of incoherently coupled NLS equations ( p^ . We assume that the number of 
negative directions (eigenvalues) of the second variation S^A is fixed, and we 
denote it as n{A). The unstable eigenvalues A of the linear problem (^) are 
connected with some negative eigenvalues of the matrix U defined by the 
elements 

The matrix U is the Hessian matrix of the soliton energy surface As{(3). We 
denote the number of positive eigenvalues of the matrix U as p{U), and the 
number of its negative eigenvalues as n{U), so that p{U) + n{U) < N, since 
some eigenvalues may be zeros in a degenerate (bifurcation) case. As is shown 
in Ref. both p{U) and n{U) satisfy some additional constraints, 

p{U) < mm{N, n{A)}, n{U) > max{0, N - n{A)}. (36) 

Within these notations, the following results on stability and instability 
of multi-component fundamental solitary waves of the coupled NLS equa- 
tions ( ^ ) can be formulated and proved |]l9| : 

• the linear problem (^2|) may have maximum n{A) unstable eigenvalues A, 
all real and positive; 

• a multi-component fundamental soliton is linearly unstable provided p(C/) < 
n{A); then the linear problem (|3^ ) has n{A) — p{U) real (positive or zero- 
becoming-positive) eigenvalues A; 

• a multi-component fundamental soliton is linearly stable provided p{U) = 
n{A){< N); in the case n{A) — N this criterion implies that the energetic 
surface As{(3) is concave up in the /3-space; 

• a single eigenvalue A crosses a marginal stability curve when the matrix U 
possesses a zero-becoming-negative eigenvalue; the normal form for the 
instability-induced dynamics of multi-component solitary waves resem- 
bles the equation of motion for an effective classical particle subjected to 
a A^-dimensional potential field, 

^ N N 

^=5i:i:«..".^^+«'(A.'). (37) 

n— 1 m— 1 

Here AInm cire the elements of the positive-definite "mass matrix" , 

N 



rfx^Gfe„(x)G'fe,„(x), 
fe=i 
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#fc(x) VJo ' df3n J 

V is the vector describing a perturbation to the soliton parameters /3, and 
W{(3, u) is an effective potential energy defined as 

N 

Wi(3, u) = AH{i3) + {f3n + i^n) APn{f3), 

n=l 

where AH {(3) =Hs{,l3 + u)- H,{(3), and Z\F„(/3) = P,„(/3 + i^) - P,„(/3), 
where index "s" marks the values calculated on the energy surface. A proof of 
those results and their comparison with the theorems of Grillakis et al. ||^,|| 
are presented in Ref. Q . 



7.2 Example: two coupled NLS equations 

In order to demonstrate how this general theory can be applied to a partic- 
ular physical problem, we consider the case of two incoherently coupled NLS 
equations, 

.dtpi 



I- 



dz dx 
. d'ip2 d'^Tp: 



dz dx^ 



(|V^2p+7lV'inV'2 = 0, (38) 



where 7 is a coupling parameter. System ( |38| ) is a two-component reduction 
of the general TV— component system ( p9| ) for di — d2 — 1, 711 = 722 = 1, and 
712 = 721 = 7- An explicit soliton solution of Eq. ( ^ ) can be easily found 
for Pi — = P and 7 > — 1 in the form, 



^^(x) =$2{x) = J-^sech(y/(3x) . (39) 
y 1 + 7 ^ ^ 

This solution describes a two-component solitary wave with the components 
of equal amplitudes, and it corresponds to a straight line (3i = P2 in the 
parameter plane (/?i,/52) of a general two-parameter family of solitary waves 
of the model (p8|). For —1 < 7 < 0, two-parameter solitons exist everywhere 
in the plane (/3i, /32), while for 7 > 0, the soliton existence region is restricted 
by two bifurcation curves, (32 — '^±i'y)Pi, where 



±2 



Approximate analytical expressions for a two-component solitary wave 
can also be obtained in the vicinity of one of the bifurcation curves (pO[), 
when one of the components of a composite solitary wave becomes small, 
while the other one is described by a scalar NLS equation. From the physical 
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point of view, this corresponds to a situation when one of the component 
creates an effective waveguide that guides the other component, and it is 
known to describe the so-caUed pulse shepherding effect in optical fibers where 
the large-ampHtude component plays a role of a shepherding pulse |29j. The 
composite soliton, that describes, for example, a large pulse V'l guiding a 
small pulse ip2, can be foimd in the form p9|. 



<?i = Ra{x) + e^R2{x) + 0(6*), = £^1(2:) + 0(e3), 
and this solution is valid in the vicinity of the bifurcation curve, 

132 = w+(7)/3i + e'w2+(7)/3i + 0{e^). 
The main terms of the asymptotic series (|4l|),(p2[) are defined as 

i?o = y/Wisech (v^x) , Si = y^sechv^ {yWix^ 

and 



(41) 
(42) 



W2+ 



!^^dx{St + 2^RoR2Sl) 



r dxSl 

and the second-order correction R2{x) is a solution of the equation, 

-5^ - 6/3i sech^ (7^2^) R2^jRqSI. 

From the existence domain of this two-component soliton, it follows that 
^^2+{'y) > 0, for < 7 < 1, and uj2+{j) < 0, for 7 > 1. At 7 = 1 (the so- 
called Manakov model) , a family of two-parameter composite solitons becomes 
degenerated: it exists on the line f3i — P2 but, generally, it is different from the 
one-parameter solution (^9|). The coupled solitons are known to be stable for 
the integrable case 7 = 1. Here we apply the stability theory developed above 
and prove that the (l-l-l)-dimensional two-parameter solitons, including the 
solitons of equal amplitudes (|39|) , ar e stable for 7 > 0, and unstable for 7 < 0. 

Following Ref. |l^ and Sec. 7.1, we evaluate the indices p{U) and n{A) 
for the explicit solution (|3^). As follows from Eqs. (|3|) and (|39|), the Hessian 
matrix U with the elements (p5f) can be found in the form. 



dPi 

dl3i 



d(32 



1 



\/^(l+7) 



and 



dPi 
dl32 



dP2 
dpi 



7 



V^(l+7)' 



Therefore, the Hessian matrix has p{U) = 2 positive eigenvalues, for —1 < 
7 < 1, and p{U) = 1 positive eigenvalue, for 7 > 1. On the other hand, the 
linear matrix operator Li given below Eq. (|l|) can be diagonalized for linear 
combinations of the eigenfunctions, vi — ui + U2 and V2 = ui — U2, such that 



(3-QI3 sech^ [^/Px^ 



Vl = ^Ui, 



«2 = ^J'V2■ 



(43) 
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Both the operators in Eqs. (^3|) are the hnear Schrodinger operators with 
solvable sech-type potentials, and the corresponding eigenvalue spectra are 
well studied. The first operator always has a single negative eigenvalue for 
fj, = —3/?, whereas the second operator has no negative eigenvalues, for 7 > 1, 
has a single negative eigenvalue, for < 7 < 1, and it has two negative 
eigenvalues, for — 1 < 7 < 0. Thus, in total there exist n{A) = 3 negative 
eigenvalues, for — 1 < 7 < 0, n{A) = 2 negative eigenvalues, for < 7 < 1, 
and n{A) — 1 negative eigenvalue, for 7 > 1. 

Applying the stability and instability results of Ref. [^| summarized in 
Sec. 7.1, we come to the conclusion that the soliton solution ( p9| ) with equal 
amplitudes is linearly stable for 7 > 0, since in this domain p{U) = n{A) = 
{1,2}, and it is linearly unstable for —1 < 7 < 0, since in this domain 
p{U) = 2 < n{A) = 3. 

Finally, we mention the other limiting case that describes the shepherding 
effect [see Eqs. ( ^ ) and (^2|)]. In this case, the elements (^) of the Hessian 
matrix U can also be calculated in an explicit analytical form, and it can 
be shown that the Hessian matrix U has p{U) = 2 positive eigenvalues, for 
< 7 < 1, and p{U) = 1 positive eigenvalue, for 7 > 1. 

On the other hand, the linear matrix operator Li cannot be diagonalized 
for the soliton in the shepherding regime (41), unless e = 0. In the latter, 
decoupled, case, it has a single negative eigenvalue at fi = — 3/3i and a double 
degenerate zero eigenvalue. When e ^ 0, the zero eigenvalue shifts to become 
H = — 2a;2+(7)/3ie^ + 0(e''). Therefore, the matrix operator Li for the soli- 
ton ( plj ) has n{A) = 2 negative eigenvalues, for < 7 < 1, and n{A) = 1 
negative eigenvalue, for 7 > 1. Thus, we come to the conclusion that the soli- 
ton in the shepherding regime is stable for 7 > since p(U) — n{A) = {1, 2}. 



8 Mult i- hump solitary waves 

In the simplest spatial optical soliton is created by one beam of a cer- 

tain polarization and frequency, and it can be viewed as a self-trapped mode 
of an effective waveguide it induces in a bulk medium Jsof . When a spatial 
soliton is composed of two (or more) modes of the induced waveguide , its 
structure becomes rather complicated, and the soliton intensity profile may 
display several peaks. This happens when one of the components creates an 
effective waveguide that guides higher-order modes of the other component. 
The resulting solitary waves are usually referred to as multi-hump solitary 
waves. Their modal structure is more complex than that of multi-component 
fundamental solitons considered in Sec. |^. 

For a long time, it was believed that all types of multi-hump solitary waves 
are linearly unstable, except for the special case of neutrally stable solitons 
in the integrable Manakov model. On the contrary, the experimental work of 
Mitchell et al. reported the observation of stationary structures resem- 



bling multi-hump solitary waves. Moreover, the recent analysis |13 confirmed 
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that multi-hump solitons supported by incoherent interaction of two optical 
beams in a photorefractive medium can indeed be stable. Below, we address 
this issue following the original results by Ostrovskaya et al. p^ . 

In the experiments , spatial multi-hump solitary waves were generated 
by incoherent interaction of two optical beams in a biased photorefractive 
medium. The corresponding model is described by a system of two coupled 
nonlinear equations for the normalized beam envelopes, u{x, z) and w{x, z), 



.du Id^u u{\u\' + \w\') 

1 -I- -4- — II — 

5Z 2 fe2 1 + s(|u|2 -f |u;|2) 

.dw Id'^w -I- jwp) 

* 97 2 9^ ^ l + s(|ii|2 + |,«|2) - = 0' 



(44) 



where the transverse, x, and propagation, z, coordinates are measured in the 
units of {Ld/ky^"^ and Ld, respectively, where Ld is a diffraction length and k 
is the wavevector. The parameter A has a meaning of a ratio of the nonlinear 
propagation constants (3i and /?2 of two components, and s is an effective 
saturation parameter (see also Sec. |l^ below). 

We look for stationary, z-independent, solutions of Eqs. ( ^ ) with both 
components u{x) and w{x) real and vanishing as \x\ — > oo. Different types of 
such two-component localized solutions, existing for < {A, s} < 1, can be 
characterized by the total power, P(A, s) = Pu +Pw, where P„ = Iwpdx 
and = \w\'^dx. If one of the components is small (e.g., w/u ^ e), 
Eqs. (^) become decoupled and, in the leading order, the equation for the 
u-component has a solution uo{x) in the form of a fundamental, sec/i-like, 
soliton with no nodes. The second equation can then be considered as an 
eigenvalue problem for the "modes" Wn{x) of a waveguide created by the 
soliton uo{x) with the effective refractive index profile ul{x)/[l + suq{x)]. Pa- 
rameter s determines the total number of guided modes and the cut-off value 
for each mode, A„(s). Therefore, a two-component vector soliton {uo,Wn) 
consists of a fundamental soliton and a n-order mode of the waveguide it 
induces in the medium, and it can be characterized by its "state vector" : 
|0,n). 

On the P(A) diagram, continuous branches representing |0,7i) solitons 
emerge at the points of bifurcations A„(s) of one-component solitons. Fami- 
lies of two-component composite solitons are found by numerical relaxation 
technique and some results of these calculations are presented in Fig. ||, 
for |0, 1) and |0, 2) solitons at s = 0.8. Observing the modification of soliton 
profiles with changing A (see inset in Fig. ^), one can see that the modal 
description of two-component solitons is valid only near the points of bifur- 
cations. For A 3> A„, the amplitude of an initially small w-component grows 
and the soliton- induced waveguide deforms. Two- and three- hump solitons 
are members of the soliton families |0, 1) (branch A-B-C) and |0,2) (branch 
D-E-F) originating at different bifurcation points. At A ^ A„(s), while the 
w-component remains small, all \0,n) solitons are single-humped, as shown 
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Fig. 6. Soliton bifurcation diagram at s = 0.8. Horizontal line - fundamental 
M— soliton. Curve A-B-C — branch of |0, 1) solitons. Curve D-E-F — branch of 
|0, 2) solitons. Inset: Transverse profiles of u (thin) and w (dashed) fields, and total 
intensity (thick), shown for marked points 




in Figs. 6(a,d). As the amplitude of the component w grows with increasing 
A, the total intensity profile, I{x) = u^ix) + w'^{x), develops [n + 1) humps 
[see Figs. 6(b,e)], and at sufficiently large A the u-component itself becomes 
multi-humped [Figs. 6(c,f)]. Separation distance between the soliton humps 
tends to infinity as A ^ 1. 

To analyze the linear stability of these multi-hump solitons, we seek solu- 
tions of Eqs. ( ^ ) in the form of weakly perturbed solitary waves: u{x, z) = 
Uq{x) -\-e[Fu{x, z) -\-iGu{x, z)] and w{x, z) = Wn{x) -|-e[F^(a:, z) + iGw{x, z)], 
where e ^ 1. Setting Fu^w ~ fu,w{x)e^^ , Gu,w ~ gu,w{x)e^^ , one can obtain 
the following eigenvalue problems 

£i£og - -Ag, (45) 
Here g = {gu,9wV, f = {fuJwV, A = y?, and 

2 (ia;2 



■^0,1 — 




1 4- X r 

; 77:72 + — Co,i 



where ag = co = //(I + si), bo = 0, ai = oq + 2mo/(1 + s/)^, ci = cq + 
2wl/{l + s/)2, and h ^ -2wow„/(l + sl)^. 

Note that CiCq and CqCi are adjoint operators with identical spectra. 
Therefore, we can consider the spectrum of only one of this operators, e.g. 
CiCq- Considering the complex yl-plane, it is straightforward to show that 



26 



Yuri S. Kivshar and Andrey A. Sukhorukov 




0.0 L I. ^ : 

0.0 0.2 0.4 0.6 0.8 1 

Soliton parameter, ^ 



Fig. 7. Existence and stability do- 
mains for two- and three-hump soli- 
tons. Shown are the existence thresh- 
olds Ai(s) and A2(s) for |0, 1) and 
|0,2) soliton families. Dashed — the 
line where total intensity of |0, 1) soli- 
ton develops humps. Shaded — insta- 
bility domain for two-hump solitons 
determined by generalized Vakhitov- 
Kolokolov criterion. Squares and cir- 
cles — instability threshold for two- 
and three-hump solitons, correspond- 
ingly, obtained by numerical solution 
of the eigenvalue problem jl^ 



A (£ (— oo,— /i^) is a continuum part of the spectrum with unbounded eigen- 
functions. Stable bounded eigenmodes of the discrete spectrum, i.e. the soli- 
ton internal modes, can have eigenvalues only inside the gap, — < Re A < 0. 
According to the general theory outlined in Sec. |^, the presence of either pos- 
itive yl or a pair of complex conjugated A^s implies soliton instability. 

Numerical solution of the eigenvalue problem ( ^ ) shows that both types 
of solitary wave solutions can be stable in a certain region of their existence 
domain, see Fig. |^. In the case of two-hump solitons, the appearance of the 
instability is related to the fact that close to the curve where the total inten- 
sity of the soliton becomes two-humped (see Fig. 0, dashed line), a pair of 
the soliton internal modes split from the continuum spectrum into the gap. 

With the aid of analytical asymptotic technique , it is possible to show 
that a perturbation mode with small but positive eigenvalue, and therefore 
linear instability of a general localized solution {u,w), appears if the deter- 
minant J{u,w} defined as 

_ djPu, P^) _Pu dPu, Pu, dPu dPu dPu, dP^ dPu 
d{/3i,(32) 2s dX 2s dX ds dX ds 9A ' 

changes its sign. Instability threshold given by the condition J = is, in 
fact, the Vakhitov-Kolokolov stability criterion |Q generalized for the case 
of two-parameter composite solitons as described in Sec. ^ However, such a 
generalization is rigorous for fundamental solitons only, whereas here we deal 
with higher-order solitary waves. As a result, in this case the determinant 
condition does not necessarily give a threshold of leading instability and, 
therefore, the presence of other instabilities (which are not associated with 
the bifurcation condition J = and can even have larger growth rates) is 
still possible. As was shown above, this is indeed true for the three-hump 
solitons, whereas the stability condition of two-hump solitons is well defined 
by the determinant condition J = 0, as verified by numerical solution of the 
linear eigenvalue problem (see Fig. |^). 
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9 Oscillatory instabilities 

When the conditions of the Vakhitov-Kolokolov criterion Q or the theorem 
of Grillakis et al. [Qjs) and their generalizations discussed in the preceding 
sections are not satisfied, for example, when there exists more than one neg- 
ative eigenvalue of the linear eigenvalue problem in a scalar case, the study 
of soliton stability should be based on the results of numerical simulations, 
and no simple stability criterion that involves the system invariants can be 
constructed. In such a case unstable eigenvalues are complex and, therefore, 
the corresponding soliton instability is called oscillatory. Oscillatory instabili- 
ties can be usually associated with resonances between two (or more) soliton 
internal modes, or resonances between the soliton internal mode and the 
edge of the continuous radiation spectrum. The former case is schematically 
presented in Fig. ||. 




Fig. 8. Schematic presentation of one 
of the scenarios of the soliton oscilla- 
tory instability that occurs via colli- 
sion of two soliton internal modes 



One of the first examples of oscillatory instabilities of solitary waves was 
found and analyzed by Barashenkov et al. for the case of optical gap 
solitons described, in the framework of the coupled-mode theory, by two 
nonlinear coupled equations for the dimensionless amplitudes of the forward- 
and backward-propagating waves, 

where a is the cross-phase-modulation parameter, and the case ct = corre- 
sponds to the exactly integrable massive Thirring model of the field theory. 

The fact that the system ( ^ ) is integrable at ct = can be employed 
to study bifurcations of new eigenvalues from the edge of the continuous spec- 
trum in the case of small but nonzero cr, similar to the case of a scalar weakly 
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perturbed cubic NLS equation (see Sec. Additional eigenvalues correspond 
to internal modes of a gap soliton, and collision of two such eigenvalues re- 
sults in a soliton oscillatory instability characterized by a pair of complex 
eigenvalues with a positive real part. 

Existence of soliton oscillatory instabilities is usually associated with 
multi-parameter solitary waves such as optical gap solitons (in the latter 
case, two independent parameters are the soliton frequency and velocity). 
Another example of the soliton oscillatory instability was recently discussed 



by Mihalache et al. |17 for the soliton propagation in birefringent optical 
fibers in the presence of walk-off, self-phase modulation, cross-phase modula- 
tion, and parametric four-wave mixing effects. As was found, the condition of 
linear marginal stability of the two-parameter solitary waves of that model is 
not necessarily given by an explicit determinant criterion, because oscillatory 
instabilities associated with complex eigenvalues were also found to exist in 
that model. 

One more interesting example of oscillatory instabilities of solitary waves 
was found for dark solitons in the generalized Ablowitz-Ladik equation and 
the discrete NLS model of the form 

+ C{4>n+l + V^n-l) + \'4>n?A^ = 0- (47) 

It was revealed that even weak inherent discreteness of a nonlinear model 
can lead to instabilities of solitary waves and, in particular, in the model ( ^7|) 
the instability of a dark soliton on a fixed background, IV'ooP = 1, appears 
for C > Ccr ~ 0.076, and it is of the oscillatory type. Additionally, it was 
found that the oscillatory instabilities may appear in a result of two different 
scenarios, and they may occur due to either a resonance between radiation 
modes and the soliton internal mode, or a resonance between two soliton 
internal modes (as shown in Fig. S). 



10 Symmetry-breaking instabilities 

All the problems discussed above involve the perturbations acting on a soliton 
in the space of the same dimension. However, a solitary wave may also become 
unstable being subjected to the action of perturbations of higher dimensions. 
We define these cases as symmetry-breaking soliton instabilities, which can be 
further classified into two categories: transverse instabilities and modulational 
instabilities. These two categories are characterized by the type of the soliton 
evolution affected either by higher spatial dimensions (such as break-up of 
a plane soliton) or temporal modulation (such as the temporal instability of 
spatial solitons and the formation of light bullets). 

A comprehensive overview of the soliton self-focusing and transverse in- 
stabilities has been recently presented by Kivshar and Pelinovsky Q, and 
we refer to that review paper for detailed discussions of a general theory of 
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the symmetry-breaking (both transverse and modulational) soliton instabil- 
ities and different examples of the instability-induced soliton dynamics, in 
the case when a plane solitary wave is subjected to perturbations of higher 
dimensions. Less studied case is the so-called azimuthal instability of radially 
symmetric solitary waves, recently observed experimentally for a number of 
physical systems. Below we present one of the examples of such symmetry- 
breaking instabilities. 

First of all, we recall that in higher dimensions there exist different types 
of radially symmetric ring-like solitary waves of higher order, even in a scalar 
NLS equation with generalized nonlinearity |]3^ . Such modes consist of a 
bright central spot surrounded by one or more rings, and all those structures 
are known to be unstable evolving into a fundamental radially symmetric 
soliton and radiation, or displaying more complicated decay-induced collapse- 
free evolution in a saturable nonlinear medium p5|| . 

More interesting types of radially symmetric localized modes are asso- 
ciated with a nonzero angular momentum of the solitary waves carrying a 
topological charge (or "spin") to, u(r, 0) = u(r)e™'^, where m — ±1, ±2, . . .. 
Such solitary waves are known as vortex solitons (see, e.g., Ref. [p6|), and 
they were found for different types of nonlinear models (see, eg. Ref. pTj). 
All types of these radially symmetric ring- like solitary waves are also unstable 
in respect to azimuthally dependent perturbations, and the decay scenario is 
determined by the value of angular momentum. In particular, due to mod- 
ulational instability the rings with the charge m — ±1 decay into pairs of 
fundamental bright solitons ("filaments") of the opposite phase, that strongly 
repel each other and always move apart along straight trajectories satisfy- 
ing the momentum conservation js^. Different filamentation scenarios were 
observed in experiments on solitons in Rb vapors js^. 

In the scalar beam propagation, symmetry- breaking instabilities of higher- 
order radially symmetric solitary waves (with or without an angular momen- 
tum) usually lead to a soliton decay or filamentation. However, for the vec- 
torial (or multi-component) nonlinear modes, the decay scenarios of ring-like 
solitary waves become more complicated and lead to novel phenomena. One 
of the recent examples of this kind is the prediction of the existence and sta- 
bility of a novel type of robust vector soliton — a dipole-mode vector soliton, 
a radially asymmetric two-component self-trapped beam in a bulk medium 
that can be created via the symmetry-breaking instability of a vortex-mode 
composite soliton earlier predicted in Ref. ||39| . 

To discuss this phenomenon, we follow the original work and consider 
two incoherently interacting beams that propagate along the direction z in 
a bulk saturable nonlinear medium. This model corresponds, in the so-called 
isotropic approximation, to the experimentally realized self-trapped beams in 
photorefractive crystals, and in the dimensionless form it can be described by 



30 



Yuri S. Kivshar and Andrey A. Sukhorukov 



a system of two coupled NLS equations for the slowly varying beam envelopes 
El and £'2, 



dz 



E 



1,2 



1 



\Ei\ 



= 0, 



(48) 



where Aj_ is the transverse Laplacian. Stationary solutions of Eq. ( p8| ) can 
be found in the form: 



El 



where Pi and P2 are two independent propagation constants. Measuring the 
transverse coordinates in the units of \/Pi, and introducing the ratio of the 
propagation constants, A (1 — /3i)/(l — /32), from Eq. (^) we obtain the 
normalized stationary equations for u and v, 



A±u -u + uf{I) ^ 0, 
A^v -\v + vf{I) = 0, 
where /(/) = /(I + sI)-\ I = 



and the parameter s 



(49) 

I- Pi plays 



the role of an effective saturation parameter, as in Eq. ([44|), Sec. 
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Fig. 9. (a) Eigenvalue of the leading unstable mode for the vortex-mode vector 
soliton (s = 0.5), and (b) typical decay scenario of the vortex- mode soliton near cut- 
off (s = 0.65, A = 0.6). Shown are intensity distributions of the u and v components 
in the {x, y) plane [M 



Radially symmetric vortex-mode vector solitons of the model (|49|) can be 
found in the form |^ : 

u{x, y) = u(r), v{x, y) = v{r)e'""^ , 
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where m is an integer topological charge. From the physical point of view, 
such self-trapped states can be regarded as composite solitons consisting of 
a self-induced waveguide created in the field u, and the first-order doughnut- 
type guided mode in the field v, the latter resembling the structure of a 
Laguerre-Gaussian (LGq) mode of a cylindrically symmetric optical waveg- 
uide. Importantly, when such a vortex-mode composite soliton propagates 
in a nonlinear medium, it undergoes a symmetry-breaking oscillatory-type 
instability, as shown in Figs. ^(a,b). However, in a sharp contrast to a vortex 
decay scenario in a scalar NLS model, two fundamental solitons created in re- 
sult of such a break-up do not move apart but, instead, they remain trapped 
in a rotating dipole-like structure, as clearly seen in Fig. ^(b). This dipole 
structure corresponds to an asymmetric vector solitary wave of a different 
type, that can be thought of as a composite self-trapped state of a soliton- 
induced waveguide that guides a dipole mode with a shape of a Hermite- 
Gaussian (HGio) laser mode. Recent studies revealed that such dipole-mode 
solitary waves are indeed extremely robust objects that survive moderate 
and large perturbations |Q. The dipole-mode composite solitons have been 
recently created experimentally in strontium barium niobate (SBN) photore- 
fractive crystals both by a phase-imprinting technique, and as a product 
of the symmetry-breaking instability of the radially symmetric vortex-mode 
solitons [|l]j4|. 

More recently, the unique robustness of the dipole-mode soliton was demon- 
strated, theoretically and experimentally, in the study of their collisions with 
scalar solitons and other dipoles P3|. 



11 Concluding remarks 

We have presented a brief overview of the key concepts and basic results of 
the soliton stability theory applied to spatial optical solitons, emphasizing the 
most recent results on the stability of multi-component solitary waves, the 
study of oscillatory instabilities of solitary waves, and symmetry-breaking 
instabilities of solitons in higher dimensions. We have demonstrated that 
the theory of soliton instabilities is closely related to the concept of soliton 
internal modes, an important feature of solitary waves of many nonintegrable 
nonlinear models. 

The simplest and most familiar class of the soliton instabilities in the 
NLS-type nonlinear models is described by the Vakhitov-Kolokolov criterion, 
dP/dP < 0, where P is the soliton power and P is its propagation constant. 
This criterion is valid for one-component fundamental solitary waves pro- 
vided the Li operator of the linear eigenvalue problem possesses exactly one 
negative eigenvalue. We have shown that a nontrivial generalization of this 
criterion to the case of multi-component solitary waves is also possible, but 
it requires more restrictive conditions. 



32 Yuri S. Kivshar and Andrey A. Sukhorukov 

When the rigorous results of the hnear stabihty theory are not available, 
a simple asymptotic theory allows to find the condition for the marginal sta- 
bility point (or, in general, a marginal stability surface) that may be further 
verified numerically. One of such cases, when the result of an asymptotic 
theory has no rigorous proof yet, even though being verified numerically, is 
the stability of two-hump solitary waves created by incoherent interaction of 
two optical beams in a photorefractive nonlinear medium. The multi-scale 
asymptotic analysis is valid for the region of the marginal stability, and it 
is a powerful tool that also allows to study weakly nonlinear effects in the 
instability-induced soliton evolution and to describe analytically different sce- 
narios of the dynamics of unstable solitons. 

In many other cases, the eigenvalues of the corresponding linear problem 
are complex, and the study of the soliton stability should be carried out 
by the numerical analysis of the corresponding linear eigenvalue problem. 
This is the case of the so-called oscillatory instabilities of solitary waves, for 
instance those known to exist for gap solitons described by a coupled-mode 
theory. Such instabilities are not predicted by simple integral criteria, and 
they appear due to a resonance between two soliton internal modes or a 
soliton internal mode and the edge of the continuous spectrum. 

Many interesting problems of the soliton stability theory and the dynam- 
ics of unstable optical solitons can be found in higher dimensions, where the 
symmetry-breaking instability of plane or radially symmetric solitary waves 
brings new features into the soliton dynamics. In particular, we have demon- 
strated that a vortex-like solitonic mode, that always breaks up into funda- 
mental bright solitons in a scalar model, may evolve, in a vector model, into 
a rotating dipole-like structure associated with a robust dipole-mode vector 
soliton. 

The study of instabilities is an important direction in the research of 
spatial optical solitons that links the theoretical concepts of the beam self- 
trapping and stationary self-guided beam propagation with the possibility of 
experimental verifications. After all, wave instabilities are probably the most 
remarkable physical phenomena that may occur in Nature, and their study 
allows researchers to reveal many important and fascinating properties of 
nonlinear physical systems. 
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